clear all; clear global; clc; close all; tstart = datetime;

load('../Set_postestimation_use.mat');

tau_A = 0.053;
tau_B = 0.038;

tau_bA = 0.053;
tau_bB = 0.041;

Par.tau_A = tau_A;
Par.tau_B = tau_B;

Par.L_A   = 1;
Par.L_B   = 1;

Par.time_ss  = 2*Par.time_ss;
Par.time_81  = Par.time1;
Par.time_95  = 20/Par.dt;
Par.time_add = (200+Par.time_81*Par.dt)/Par.dt-Par.time_ss+1;

Par.num_r = Par.num_r/3;

[Par,Step,Trans] = Trade_ctrf(Par);

horizons = 5:5:50;
tau_bA_grid = linspace(0,0.60,301);
trfA_grid   = 0.1:-0.005:00;

rev_b   = 1; % Revenue back to consumer
retal_b = 0; % 1 if foreign retaliation

Par.rev_US = rev_b;
Par.rev_FN = rev_b;

% ============================================
%                 SOLUTIONS
% ============================================

Param = structvars(Par);
for ctr0 = 1:size(Param,1); eval(Param(ctr0,:)); end

steps      = Step{1};
Trans_innov = Step{2};
Trans_fail = Step{3};
Trans_exo  = Step{4};

%% Calibrated Economy Path
% =================================================================

Tol_rr = 1*10^-6;
save Tol_rr Tol_rr

%---------------- Path

mu_init = y_init;
Q_init  = y_init;

Par.if_long = 0;

[rr_sqn, X_sqn, V_sqn, Xe_sqn, Q_sqn, Ql_sqn, mu_share, ind_neg, ...
                            prft_fnc_sqn, wage_dom_sqn, mcuts] = ...
                                Path_ctrf(mu_init,Q_init,Par,Step,Trans);
                            
%---------------- Variables

% Profits
% --------------------------------------------------------

prft_fnc_A = prft_fnc_sqn(1:size(prft_fnc_sqn,1)/2,:);
prft_fnc_B = prft_fnc_sqn(size(prft_fnc_sqn,1)/2+1:end,:);

wage_dom_A = wage_dom_sqn(1:size(wage_dom_sqn,1)/2,:);
wage_dom_B = wage_dom_sqn(size(wage_dom_sqn,1)/2+1:end,:);

% Qualities
% --------------------------------------------------------

QA = Q_sqn(1:size(Q_sqn,1)/2,:);
QB = Q_sqn(size(Q_sqn,1)/2+1:end,:);

QA_long = Ql_sqn(1:size(Q_sqn,1)/2,:);
QB_long = Ql_sqn(size(Q_sqn,1)/2+1:end,:);

QwA = sum(wage_dom_A.*QA+flip(1-wage_dom_B).*QA*(1+kapa)/((1+kapa)*(1+trf_B))^(1/bet));
QwB = sum(wage_dom_B.*QB+flip(1-wage_dom_A).*QB*(1+kapa)/((1+kapa)*(1+trf_A))^(1/bet));

qq_sqn = [sum(QA)./QwA; sum(QB)./QwB];
                
wq_sqn = [(1-Par.bet)*Par.eta^(Par.bet-1)*qq_sqn(1,:).^(-Par.bet);
            (1-Par.bet)*Par.eta^(Par.bet-1)*qq_sqn(2,:).^(-Par.bet)];  

xk_sqn = [(Par.pi_frn/Par.bet)^(1/(1-Par.bet))./wq_sqn(1,:).^(1/Par.bet)*Par.L_B;
         (Par.pi_frn_toUS/Par.bet)^(1/(1-Par.bet))./wq_sqn(2,:).^(1/Par.bet)*Par.L_A];

% Shares
% --------------------------------------------------------

MU      = mu_share(:,1:time_ss+1);
MU_long = mu_share;

McutX  = mcuts(1,1:time_ss);
McutM  = mcuts(2,1:time_ss);

% R&D
% --------------------------------------------------------

X_iA = X_sqn(1:size(X_sqn,1)/2,:);
X_iB = X_sqn(size(X_sqn,1)/2+1:end,:);

V_iA = V_sqn(1:size(V_sqn,1)/2,:);
V_iB = V_sqn(size(V_sqn,1)/2+1:end,:);

X_eA = Xe_sqn(1:size(Xe_sqn,1)/2,:);
X_eB = Xe_sqn(size(Xe_sqn,1)/2+1:end,:);

V_eA = alfa_eA*X_eA.^gama_eA*(1-1/gama_eA);
V_eB = alfa_eB*X_eB.^gama_eB*(1-1/gama_eB);

% Income
% --------------------------------------------------------

II.IA_1_pr  = sum(prft_fnc_A.*QA).*wq_sqn(1,:).^(1-1/bet);
II.IA_2_rdi = sum(alfa_A/gama_A*X_iA.^gama_A.*QA);
II.IA_3_rde = sum(alfa_eA/gama_eA*X_eA.^gama_eA.*QA);
II.IA_4_wd  = pi_dom/(1-bet)*L_A*wq_sqn(1,:).^(1-1/bet).*sum(wage_dom_A.*QA);
II.IA_5_wf  = pi_frn_toUS/(1-bet)*L_A*wq_sqn(2,:).^(1-1/bet).*sum(flip(1-wage_dom_A).*QB);
II.IA_6_wig = wq_sqn(1,:).*sum(QA);
II.IA_7_trf = rev_US*trf_A*((1+trf_A)*(1+kapa)*eta/(1-bet))*wq_sqn(2,:).*xk_sqn(2,:).*sum(flip(1-wage_dom_A).*QB);

II.IB_1_pr  = sum(prft_fnc_B.*QB).*wq_sqn(2,:).^(1-1/bet);
II.IB_2_rdi = sum(alfa_B/gama_B*X_iB.^gama_B.*QB);
II.IB_3_rde = sum(alfa_eB/gama_eB*X_eB.^gama_eB.*QB);
II.IB_4_wd  = pi_dom/(1-bet)*L_B*wq_sqn(2,:).^(1-1/bet).*sum(wage_dom_B.*QB);
II.IB_5_wf  = pi_frn/(1-bet)*L_B*wq_sqn(1,:).^(1-1/bet).*sum(flip(1-wage_dom_B).*QA);
II.IB_6_wig = wq_sqn(2,:).*sum(QB);
II.IB_7_trf = rev_FN*trf_B*((1+trf_B)*(1+kapa)*eta/(1-bet))*wq_sqn(1,:).*xk_sqn(1,:).*sum(flip(1-wage_dom_B).*QA);

IA = sum(prft_fnc_A.*QA).*wq_sqn(1,:).^(1-1/bet)-...
    sum(alfa_A/gama_A*X_iA.^gama_A.*QA)-...
    sum(alfa_eA/gama_eA*X_eA.^gama_eA.*QA)+...
    pi_dom/(1-bet)*L_A*wq_sqn(1,:).^(1-1/bet).*sum(wage_dom_A.*QA)+...
    pi_frn_toUS/(1-bet)*L_A*wq_sqn(2,:).^(1-1/bet).*sum(flip(1-wage_dom_A).*QB)+...
    wq_sqn(1,:).*sum(QA)+...
    rev_US*trf_A*((1+trf_A)*(1+kapa)*eta/(1-bet))*wq_sqn(2,:).*xk_sqn(2,:).*sum(flip(1-wage_dom_A).*QB);  
IB = sum(prft_fnc_B.*QB).*wq_sqn(2,:).^(1-1/bet)-...
    sum(alfa_B/gama_B*X_iB.^gama_B.*QB)-...
    sum(alfa_eB/gama_eB*X_eB.^gama_eB.*QB)+...
    pi_dom/(1-bet)*L_B*wq_sqn(2,:).^(1-1/bet).*sum(wage_dom_B.*QB)+...
    pi_frn/(1-bet)*L_B*wq_sqn(1,:).^(1-1/bet).*sum(flip(1-wage_dom_B).*QA)+...
    wq_sqn(2,:).*sum(QB)+...
    rev_FN*trf_B*((1+trf_B)*(1+kapa)*eta/(1-bet))*wq_sqn(1,:).*xk_sqn(1,:).*sum(flip(1-wage_dom_B).*QA); 

g_IA = mean((IA(2:time_81+1)-IA(1:time_81))./IA(1:time_81));
g_IB = mean((IB(2:time_81+1)-IB(1:time_81))./IB(1:time_81)); 

gI_A = (IA(2:end)./IA(1:end-1))-1;
gI_B = (IB(2:end)./IB(1:end-1))-1;


% Welfare
% --------------------------------------------------------

if psy == 1
    
    for ctr_t = 1:length(horizons)
    
        Wlf_o(ctr_t,1) = sum(exp(-ro*(dt:dt:horizons(ctr_t))).*log(IA(time_81+1:time_81+horizons(ctr_t)/dt)));

        Wlf_oB(ctr_t,1) = sum(exp(-ro*(dt:dt:horizons(ctr_t))).*log(IB(time_81+1:time_81+horizons(ctr_t)/dt)));
        
    end  
    
else
    
    for ctr_t = 1:length(horizons)
        
        Wlf_o(ctr_t,1) = (1-psy)^-1*[sum(exp(-ro*(dt:dt:horizons(ctr_t))).*(IA(time_81+1:time_81+horizons(ctr_t)/dt).^(1-psy)))];

        Wlf_oB(ctr_t,1) = (1-psy)^-1*[sum(exp(-ro*(dt:dt:horizons(ctr_t))).*(IB(time_81+1:time_81+horizons(ctr_t)/dt).^(1-psy)))];
        
    end  
   
end
  

%% Policy Change
% =================================================================
%

% ------------- Parfor intialization ------------- 

mub_init = mu_share(:,time_81+1);
Qb_init  = QA(:,time_81+1);

mu_pre = mu_share(:,1:time_81);
rr_pre = rr_sqn(:,1:time_81);
IA_pre = IA(:,1:time_81);
IB_pre = IB(:,1:time_81);
QA_pre = QA(:,1:time_81);
QB_pre = QB(:,1:time_81);
QwA_pre = QwA(:,1:time_81);
QwB_pre = QwB(:,1:time_81);
McutM_pre = McutM(1,1:time_81);
McutX_pre = McutX(1,1:time_81);

QA_post = QA(:,time_81+1:end);
QB_post = QB(:,time_81+1:end);

II_orig = II;

rr_ss_run = rr_sqn(:,end)*ones(1,length(trfA_grid));
% save rr_guess_run rr_ss_run

Tol_rr = 1*10^-6;
save Tol_rr Tol_rr

% parpool(24)

for ctr_k = 1:length(trfA_grid)

    trfA_call = trfA_grid(ctr_k);

    Par_b = Par;
    % Par_b.num_r = Par.num_r/2;

    Par_b.trf_A = trfA_grid(ctr_k); 
    Par_b.trf_B = Par.trf_B+retal_b*(trfA_grid(ctr_k)-Par.trf_A);
    
    Par_b.ctr_k = ctr_k;

    [Par_b,Step_b,Trans_b] = Trade_ctrf(Par_b);

    % ============================================
    %                 SOLUTIONS
    % ============================================

    pi_domb = Par_b.pi_dom;
    pi_frnb = Par_b.pi_frn;
    pi_frnb_toUS = Par_b.pi_frn_toUS;

    Par_loop = Par_b;

    parfor ctr_s = 1:length(tau_bA_grid)
        
        Run = [ctr_s,ctr_k];

        tau_bA = tau_bA_grid(ctr_s);

        Par_b = Par_loop;
    
        Par_b.tau_A = tau_bA;
        Par_b.tau_B = tau_bB;

        %---------------- Path

        [rrb_sqn, Xb_sqn, Vb_sqn, Xeb_sqn, Qb_sqn, Qlb_sqn, mub_share, ind_neg_b, ...
                                    prft_fncb_sqn, wage_domb_sqn, mcutsb] = ...
                                    Path_ctrf(mub_init,Qb_init,Par_b,Step_b,Trans_b);

        rr_b = [rr_pre rrb_sqn(:,1:end-time_81)];
        
        rr_ss_run(:,ctr_s) = rr_b(:,end);

        %---------------- Variables

        % Profits
        % --------------------------------------------------------

        prft_fncb_A = prft_fncb_sqn(1:size(prft_fncb_sqn,1)/2,:);
        prft_fncb_B = prft_fncb_sqn(size(prft_fncb_sqn,1)/2+1:end,:);

        wage_domb_A = wage_domb_sqn(1:size(wage_domb_sqn,1)/2,:);
        wage_domb_B = wage_domb_sqn(size(wage_domb_sqn,1)/2+1:end,:);

        % Qualities
        % --------------------------------------------------------

        QA_b = Qb_sqn(1:size(Q_sqn,1)/2,:);
        QB_b = Qb_sqn(size(Q_sqn,1)/2+1:end,:);

        QAb_long = Qlb_sqn(1:size(Q_sqn,1)/2,:);
        QBb_long = Qlb_sqn(size(Q_sqn,1)/2+1:end,:);

        QwA_b = sum(wage_domb_A.*QA_b+flip(1-wage_domb_B).*QA_b*(1+kapa)/((1+kapa)*(1+Par_b.trf_B))^(1/bet));
        QwB_b = sum(wage_domb_B.*QB_b+flip(1-wage_domb_A).*QB_b*(1+kapa)/((1+kapa)*(1+Par_b.trf_A))^(1/bet));

        qqb_sqn = [sum(QA_b)./QwA_b; sum(QB_b)./QwB_b];

        wqb_sqn = [(1-Par_b.bet)*Par_b.eta^(Par_b.bet-1)*qqb_sqn(1,:).^(-Par_b.bet);
                    (1-Par_b.bet)*Par_b.eta^(Par_b.bet-1)*qqb_sqn(2,:).^(-Par_b.bet)];  

        xkb_sqn = [(Par_b.pi_frn/Par_b.bet)^(1/(1-Par_b.bet))./wqb_sqn(1,:).^(1/Par_b.bet)*Par.L_B;
                    (Par_b.pi_frn_toUS/Par_b.bet)^(1/(1-Par_b.bet))./wqb_sqn(2,:).^(1/Par_b.bet)*Par.L_A];

        QA_orgb = [QA_pre QA_b(:,1:end-time_81)];
        QB_orgb = [QB_pre QB_b(:,1:end-time_81)];

        QwA_orgb = [QwA_pre QwA_b(:,1:end-time_81)];
        QwB_orgb = [QwB_pre QwB_b(:,1:end-time_81)];

        qq_orgb  = [sum(QA_orgb)./QwA_orgb; sum(QB_orgb)./QwB_orgb];

        % Shares
        % --------------------------------------------------------

        MU_b = [mu_pre mub_share(:,1:end-time_81)];

        MCUTX(ctr_k,ctr_s).cut  = [McutX_pre mcutsb(1,1:time_ss-time_81)];
        MCUTM(ctr_k,ctr_s).cut  = [McutM_pre mcutsb(2,1:time_ss-time_81)];

        % R&D
        % --------------------------------------------------------

        Xb_iA = Xb_sqn(1:size(X_sqn,1)/2,:);
        Xb_iB = Xb_sqn(size(X_sqn,1)/2+1:end,:);

        Xb_eA = Xeb_sqn(1:size(Xe_sqn,1)/2,:);
        Xb_eB = Xeb_sqn(size(Xe_sqn,1)/2+1:end,:);

        % Income 
        % --------------------------------------------------------

        II = II_orig;

        II.IA_bb1_pr  = sum(prft_fncb_A.*QA_b).*wqb_sqn(1,:).^(1-1/bet);
        II.IA_bb2_rdi = sum(alfa_A/gama_A*Xb_iA.^gama_A.*QA_b);
        II.IA_bb3_rde = sum(alfa_eA/gama_eA*Xb_eA.^gama_eA.*QA_b);
        II.IA_bb4_wd  = pi_domb/(1-bet)*L_A*wqb_sqn(1,:).^(1-1/bet).*sum(wage_domb_A.*QA_b);
        II.IA_bb5_wf  = pi_frnb_toUS/(1-bet)*L_A*wqb_sqn(2,:).^(1-1/bet).*sum(flip(1-wage_domb_A).*QB_b);
        II.IA_bb6_wig = wqb_sqn(1,:).*sum(QA_b);
        II.IA_bb7_trf = Par_b.rev_US*Par_b.trf_A*((1+Par_b.trf_A)*(1+kapa)*eta/(1-bet))*wqb_sqn(2,:).*xkb_sqn(2,:).*sum(flip(1-wage_domb_A).*QB_b);

        II.IA_bb2_xi = sum(alfa_A/gama_A*X_iA.^gama_A.*[QA(:,1:time_81) QA_b(:,1:end-time_81)]);
        II.IA_bb3_xe = sum(alfa_eA/gama_eA*X_eA.^gama_eA.*[QA(:,1:time_81) QA_b(:,1:end-time_81)]);

        II.IB_bb1_pr  = sum(prft_fncb_B.*QB_b).*wqb_sqn(2,:).^(1-1/bet);
        II.IB_bb2_rdi = sum(alfa_B/gama_B*Xb_iB.^gama_B.*QB_b);
        II.IB_bb3_rde = sum(alfa_eB/gama_eB*Xb_eB.^gama_eB.*QB_b);
        II.IB_bb4_wd  = pi_domb/(1-bet)*L_B*wqb_sqn(2,:).^(1-1/bet).*sum(wage_domb_B.*QB_b);
        II.IB_bb5_wf  = pi_frnb/(1-bet)*L_B*wqb_sqn(1,:).^(1-1/bet).*sum(flip(1-wage_domb_B).*QA_b);
        II.IB_bb6_wig = wqb_sqn(2,:).*sum(QB_b);
        II.IB_bb7_trf = Par_b.rev_FN*Par_b.trf_B*((1+Par_b.trf_B)*(1+kapa)*eta/(1-bet))*wqb_sqn(1,:).*xkb_sqn(1,:).*sum(flip(1-wage_domb_B).*QA_b);

        II.IB_bb2_xi = sum(alfa_B/gama_B*X_iB.^gama_B.*[QB(:,1:time_81) QB_b(:,1:end-time_81)]);
        II.IB_bb3_xe = sum(alfa_eB/gama_eB*X_eB.^gama_eB.*[QB(:,1:time_81) QB_b(:,1:end-time_81)]);

        IA_bb = II.IA_bb1_pr-II.IA_bb2_rdi-II.IA_bb3_rde+...
                II.IA_bb4_wd+II.IA_bb5_wf+II.IA_bb6_wig+II.IA_bb7_trf;

        IB_bb = II.IB_bb1_pr-II.IB_bb2_rdi-II.IB_bb3_rde+...
                II.IB_bb4_wd+II.IB_bb5_wf+II.IB_bb6_wig+II.IB_bb7_trf;

        IA_b = [IA_pre IA_bb(:,1:end-time_81)];
        IB_b = [IB_pre IB_bb(:,1:end-time_81)];

        % Alternative for Welfare

        QA_bo = [QA_post 0*QA_b(:,end-time_81+1:end)];
        QB_bo = [QB_post 0*QB_b(:,end-time_81+1:end)];

        QwA_bo = sum(wage_domb_A.*QA_bo+flip(1-wage_domb_B).*QA_bo*(1+kapa)/((1+kapa)*(1+Par_b.trf_B))^(1/bet));
        QwB_bo = sum(wage_domb_B.*QB_bo+flip(1-wage_domb_A).*QB_bo*(1+kapa)/((1+kapa)*(1+Par_b.trf_A))^(1/bet));

        qqbo_sqn = [sum(QA_bo)./QwA_bo; sum(QB_bo)./QwB_bo];

        wqbo_sqn = [(1-Par.bet)*Par.eta^(Par.bet-1)*qqbo_sqn(1,:).^(-Par.bet);
                    (1-Par.bet)*Par.eta^(Par.bet-1)*qqbo_sqn(2,:).^(-Par.bet)];  

        xkbo_sqn = [(Par.pi_frn/Par.bet)^(1/(1-Par.bet))./wqbo_sqn(1,:).^(1/Par.bet)*Par.L_B;
                    (Par.pi_frn_toUS/Par.bet)^(1/(1-Par.bet))./wqbo_sqn(2,:).^(1/Par.bet)*Par.L_A];

        II.IA_bo1_pr  = sum(prft_fncb_A.*QA_bo).*wqbo_sqn(1,:).^(1-1/bet);
        II.IA_bo2_rdi = sum(alfa_A/gama_A*Xb_iA.^gama_A.*QA_bo);
        II.IA_bo3_rde = sum(alfa_eA/gama_eA*Xb_eA.^gama_eA.*QA_bo);
        II.IA_bo4_wd  = pi_domb/(1-bet)*L_A*wqbo_sqn(1,:).^(1-1/bet).*sum(wage_domb_A.*QA_bo);
        II.IA_bo5_wf  = pi_frnb_toUS/(1-bet)*L_A*wqbo_sqn(2,:).^(1-1/bet).*sum(flip(1-wage_domb_A).*QB_bo);
        II.IA_bo6_wig = wqbo_sqn(1,:).*sum(QA_bo);
        II.IA_bo7_trf = Par_b.rev_US*Par_b.trf_A*((1+Par_b.trf_A)*(1+kapa)*eta/(1-bet))*wqbo_sqn(2,:).*xkbo_sqn(2,:).*sum(flip(1-wage_domb_A).*QB_bo);

        II.IB_bo1_pr  = sum(prft_fncb_B.*QB_bo).*wqbo_sqn(2,:).^(1-1/bet);
        II.IB_bo2_rdi = sum(alfa_B/gama_B*Xb_iB.^gama_B.*QB_bo);
        II.IB_bo3_rde = sum(alfa_eB/gama_eB*Xb_eB.^gama_eB.*QB_bo);
        II.IB_bo4_wd  = pi_domb/(1-bet)*L_B*wqbo_sqn(2,:).^(1-1/bet).*sum(wage_domb_B.*QB_bo);
        II.IB_bo5_wf  = pi_frnb/(1-bet)*L_B*wqbo_sqn(1,:).^(1-1/bet).*sum(flip(1-wage_domb_B).*QA_bo);
        II.IB_bo6_wig = wqbo_sqn(2,:).*sum(QB_bo);
        II.IB_bo7_trf = Par_b.rev_FN*Par_b.trf_B*((1+Par_b.trf_B)*(1+kapa)*eta/(1-bet))*wqbo_sqn(1,:).*xkbo_sqn(1,:).*sum(flip(1-wage_domb_B).*QA_bo);

        % Welfare
        % --------------------------------------------------------

        [WlfA,WlfB,IAb_set,IBb_set] = postestimation_welfare_trade(II,Par_b,horizons); 

        WLFA(:,ctr_s) = WlfA(:,1);
        WLFB(:,ctr_s) = WlfB(:,1);

        WLF_all(ctr_k,ctr_s).rateA = [trfA_call tau_bA_grid(ctr_s)];
        WLF_all(ctr_k,ctr_s).subsA = WlfA;
        WLF_all(ctr_k,ctr_s).subsB = WlfB;
        WLF_all(ctr_k,ctr_s).flag  = ind_neg_b;

    end    

    CEQ(:,ctr_k,:) = (WLFA./(Wlf_o*ones(1,size(WLFA,2)))).^(1/(1-psy));  

    WLFA_both(:,ctr_k,:) = WLFA;
    WLFB_both(:,ctr_k,:) = WLFB;
    
    % save rr_guess_run rr_ss_run

end

% save Optimal_combined_over_time4 CEQ tau_bA_grid trfA_grid horizons
% save Optimal_combined_over_time_data4 CEQ WLF_all WLFA_both WLFB_both Wlf_o psy tau_bA_grid trfA_grid horizons McutX McutM MCUTX MCUTM
% 
% save Optimal_rnd_over_open4.mat CEQ tau_bA_grid trfA_grid horizons
% save Optimal_rnd_over_open_data4.mat CEQ WLF_all WLFA_both WLFB_both Wlf_o psy tau_bA_grid trfA_grid horizons McutX McutM MCUTX MCUTM

% save Open_deneme_init1 CEQ tau_bA_grid trfA_grid horizons

dur = datetime-tstart; fprintf('This message is sent at time %s ending the job that started at time %s and took %s long.\n', datetime, tstart, dur);

delete(gcp('nocreate'))
